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Abstract 

By the standard theory, the stable Qk+i,k-Qk,k+i/Qi c divergence-free element con- 
verges with the optimal order of approximation for the Stokes equations, but only order k 
for the velocity in iT -norm and the pressure in L 2 -norm. This is due to one polynomial 
degree less in y direction for the first component of velocity, a Qk+i,k polynomial. In this 
manuscript, we will show a superconvergence of the divergence free element that the order 
of convergence is truly k + 1, for both velocity and pressure. Numerical tests are provided 
confirming the sharpness of the theory. 
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1 Introduction 

The divergence- free finite element method is mainly for solving incompressible flow problems, 
such as Stokes or Navier-Stokes equations, where the finite element space for the pressure is 
exactly the divergence of the finite element space for the velocity. In such a method, the finite 
element velocity is divergence-free pointwise, i.e. the incompressibility condition is enforced 
strongly. Traditional finite elements enforce the incompressibility weakly, cf. [17\ [9]. That is, 
in order to satisfy the inf-sup stability condition, the incompressibility condition is weakened 
by either enlarging the velocity space or decreasing the pressure space. This often leads to 
some sub-optimal methods, or a waste of computation, due to the imperfect matching of two 
spaces. It may lead to inaccurate mass conservation, which is critical in certain computational 
problems. 

A fundamental study on the divergence- free element method was done by Scott and Vogelius 
( \18\ I19j ) that the Pfc+i/P/r method is stable and consequently of the optimal order on 2D 
triangular grids, for k > 3. Here the velocity space is the continuous piecewise-polynomials of 
degree (k + 1) or less while the the pressure space is the discontinuous piecewise-polynomials 
of degree k or less, or the divergence of the velocity, to be precise. There are several other such 
divergence-free finite elements, cf. [21 QH H5J QH [§3j [Ml I2S] • 

Starting from the most popular element, the Qi/Pq element ([BJ [7j), there is a series of 
work on mixed finite elements on rectangular grids in 2D and 3D. Brezzi and Falk showed 
that the Qk+i/Qf c element is unstable in [10], for any k > 0. Here Qf c denotes the space of 
discontinuous piecewise-polynomials. In [21], Stenberg and Suri showed the stability, but a sub- 
optimal order of approximation, for the Qk+i/QkL\ element for all k > 1 in 2D. Bernardi and 
Maday proved the stability and the optimal order of convergence for the Qk+\/P^ c element, cf. 
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[3] . Ainsworth and Coggins established pQ the stability and the optimal order of convergence for 
the Taylor-Hood Qk+i/Qk element, where the pressure space is continuous too. The Bernardi- 
Raugel element ([5]) optimizes the Qk+i/Q^Li element, when k = 1, by reducing the velocity 
space to Qi 2-Q2 1 polynomials. Here the first component of velocity in the Bernardi-Raugel 
element is a polynomial of degree 1 in x direction, but of degree 2 in y direction. To be precise, 
the Bernardi-Raugel element enrich the Q\ velocity space by face-bubble functions. Similar 
to the Bernardi-Raugel element, a divergence-free finite element, Qk+i,k-Qk,k+i/Qk C (^ > 2), 
was proposed in [25j . which further optimizes the Bernardi-Raugel element by increasing the 
polynomial degree of pressure from (A;— 1) to k. The nodal degrees of freedom of this divergence- 
free element and the Bernardi-Raugel element are plotted in Figure [TJ This divergence-free 
element was extended to its lowest-order form, k = 1, i.e., Q2,i-Qi,2/Qf c , in [II]- Here the 
space Qf for the pressure is the space of discontinuous Qk polynomials with all spurious 
modes removed, i.e., eliminating one degree of freedom at each vertex. In the construction, 
the pressure space is exactly the divergence of the velocity. Thus, the resulting finite element 
is divergence-free pointwise. In such a case, the discrete pressure space can be omitted in the 
computation. By an iterated penalty method, we obtain the pressure solution as a byproduct, 
cf. [24] and Section |4] below. However, by the standard finite element theory developed in 
[141 125] , this divergence-free element converges at order k only, due to a degree k polynomial 
in y for the first component of u/j. This cannot be improved by the standard theory, where 
the optimal order of convergences is derived from the inf-sup stability. In this manuscript, we 
further study this Qk+ i,fe-Qfe,fc+idivergence-free element and show its superconvergence, that it 
does converge at order k + 1. Further the velocity of the Qfc+i^-Qfc^+idivergence-free element 
may be ultraconvergent, i.e., two orders higher than the standard convergence, provided the 
interpolation polynomial is divergence-free. The extension of this divergence-free element to 
3D is straightforward, so is its superconvergence property. 




Ph- 



Ph- 




Figure 1: Nodes of u^/ph for divergence-free (top) and Bernardi-Raugel elements. 



The rest of the paper is organized as follows. In Section 2, we define the finite element 
for the Stationary Stokes equations. In Section 3, we establish a superconvergence for the 
divergence-free element. In Section HI we provide some test results confirming the analysis. In 
particular, we show the order of convergence of the divergence-free element is one higher than 
that of the rotated Bernardi-Raugel element. 



2 



2 The Qk+i,k~Qk,k+i divergence-free element 

In this section, we shall define the divergence-free finite element for the stationary Stokes 
equations on rectangular grids. The resulting finite element solutions for the velocity are 
divergence-free point wise. 

We consider a model stationary Stokes problem: Find the velocity u and the pressure p on 
a 2D polygonal domain $7, which can be subdivided into rectangles, such that 

-Au + Vp = f in O, 

divu = infi, (2.1) 

u = on dfl. 

The weak form for flUJ is: Find u G H^{Q) 2 and p G L§(fi) := L 2 {Q)/C = {p G L 2 | J n p = 0} 

a(u,v)+6(v,p) = (f,v) VvG# W, 
6(u,g) = \/q G £,§(«). 



Here i?o(J7) 2 is the subspace of the Sobolev space ff 1 (0) 2 (cf. [TT]) with zero boundary trace, 
and 



a(u, v) = / Vu • Vv dx, 
isi 

^( V )P) = — / divvpefcc, 
(f , v) = / f v dx. 



Figure 2: Three levels of grids, and a macro-element grid (for k = 1 only). 

The finite element grids are defined by, cf. Figure 

% = {K | UK = J2, K = [x a , x b ] x [y c , y d ] with size /i^ = m&x{x b - x a , y d - y c } < h] . 

We further assume, only for the lowest-order element k = 1 in (12.3H . that the rectangles in grid 
7/,, can be combined into groups of four to form a macro-element grid: 

M h = {M\M = uf =1 Ki = [xi-i,x i+1 ] x [ife-i.jfc+i], if, G T h , = 0} . 

See the 4th diagram in Figure [2) The polynomial spaces are defined by 

Qk,i = { cy^y'}, Q fc = 

i<k,j<l 
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The Qk+i,k-Qk,k+i(k > 1) element spaces are 

V ft = {v /t € C(n) 2 | v h \ K G Qk+i,k x Qk,k+\ G T h , and u^| an = 0} , (2.3) 
P ft = {divu fo | u ft G V A } . (2.4) 

Since J* n p/i = J n div = J dn = for any ph G Ph, we conclude that 

v h cH^nf, p h cL 2 (n), 

i.e., the mixed- finite element pair is conforming. The resulting system of finite element equa- 
tions for (|2.2|) is: Find G V/j and ph G -f\ such that 

a(u h ,v) + b(v,p h ) = (f,v) VvGV fc , 
6(u h ,g) = VgeP*. 

Traditional mixed-finite elements require the inf-sup condition to guarantee the existence 
of discrete solutions. As (j2.4j) provides a compatibility between the discrete velocity and the 
discrete pressure spaces, the linear system of equations (|2.5p always has a unique solution, cf. 
[24] . Furthermore, such a solution u/j is divergence-free: by the second equation in (|2.5p and 
the definition of P h in (l2T4"jl . 

&( u ft><?) = K u /i>- divu /i) = II div ^11^(^)2 = 0. (2.6) 

In this case, i.e., V/j C Z := {divv | v G Hq(Q,) 2 }, we call the mixed finite element a 
divergence- free element. It is apparent that the discrete velocity solution is divergence- free if 
and only if the discrete pressure finite element space is the divergence of the discrete velocity 
finite element space, i.e., (12. 4j) . 

We note that by (|2.4p . -f\ is a subspace of discontinuous, piecewise bilinear polynomials. 
As singular vertices are present (see [TBI EHJ El ES] ) , -P/i is a proper subset of the discontinuous 
piecewise Q\ polynomials. It is possible, but difficult to find a local basis for P^. But on 
the other side, it is the special interest of the divergence-free finite element method that the 
space Ph can be omitted in computation and the discrete solutions approximating the pressure 
function in the Stokes equations can be obtained as byproducts, if an iterated penalty method 
is adopted to solve the system (|2.5p . cf. [131 El El 1201 [2l] for more information. 

3 Superconvergence 

As usual, the superconvergence is obtained by the method of integration by parts, cf. [l~2l 122]. 
But we have a long series of lemmas dealing with each term in the bilinear forms a(-, ■) and 

For a convenience in referring components of the vector velocity, we define the two inho- 
mogeneous polynomial spaces: 

V h ,i = {cp G H^(n) | <j)\ K G Q k+ i,k Vif G T h }, (3.1) 
V h ,2 = W G H^) | 4>\k G Q k ,k+i VK G T h }, (3.2) 

k > 1. That is, 

V h = V K i x V h ,2, k>l. 
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® u d a i ) = u ( a i ) 

. J uip 2 ds = J up 2 ds 

o / uiPids = J upids 

m J K uiq 2 ,idx = J K uq 2 ,idx 

o j K uiqi t2 dx = J K uq lt2 dx 



Q?,,2 Q2,3 

Figure 3: Three types of interpolation nodes, k = 3. 



The interpolation operator \% is defined for the two components of u: 

I h :H^(n)xH^)^V hA xV hi2 , 



(3.3) 



To define uj at the Lagrange nodes, we define its vertex nodal values, then internal edge values, 
and finally internal values, by solving the following equations sequentially (see Figure [3]): 



(u-ui){af) 


= 


at four vertices of K, MK € Th, 


(3.4) 


\ (u - ui)p k _i(x)dx 


= 


on the top and bottom edges of K, 


(3.5) 


h=vj 








J (u - ui)p k -2{y)dy 


= 


on the left and right edges of K. 


(3.6) 


/ (u - ui)q k -i,k-2dx 

JK 


= 


on the square K, 


(3.7) 



where p^ G the space of ID polynomials of degree k or less, and q^j € Qk,l- By rotating x 
and y, vi is defined similarly /symmetrically to ui. 



Lemma 3.1 (two-order superconvergence) For any Qk+ik function if) € V/j \, defined in (|3.1 
for any u G H k+3 (£l), and for all k > I, 



I (u- ui) x ip x dx\ < Ch k+2 \\u\\ Hk +3\\ip\\ H i. 
Jn 



(3.8) 



Proof. We first consider the estimation on the reference element K = [— 1,1] . Since ij) G 
Qk+i jfe, we have an exact Taylor expansion: 



yk 1 yk 

ip x (x,y) = ip x (x,0) + yi/j xy (x,0) H h ^ _ ^ i() xy k-i(x,0) + — ip xy k(x, 0), 



(3.9) 



where ip x (x, 0) and all ip xy j (x, 0) are Pk polynomials in x only. We will perform the integration 
by parts repeatedly. First, for the lower order terms in (|3.9p . we notice that, by the definition 
of ui in ([33]) and (1377) . 



/ (n - ui) x y J ip xy j (x, 0)dx 

JK 



(u - uijy il) xyj {x,ti)\ x x= _ x dy - l(u- u I )y J ip x2yJ (x,0)dx 



K 



=0 



when j = 0, 1, . . . , k — 2. 



(3.10) 
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Please be aware that tp x 2 y j(x,0) E Pk-\(x). Hence, we need to deal with only the last two 
terms in (I3.9p . 

For the last two terms in (j3.9|) . in order to do integration by parts, we express polynomials 
y k ~ 1 and y k by derivatives of another polynomial. 



Sk(y) 



( y 2 _ 1)fc +i y 2k+2 



2k 



(2k + 2)\ {2k + 2)\ (2k + 2)\ 
4 J) (±1) = 0, j =0,1,-.. , A;, 



+ 



y 



2k+2 



(2k + 2)1 



s[ k+2 \y) = ^y k + P k-2(y), 



+ P2 k (y), (3.11) 

(3.12) 
(3.13) 



Here p2k(y) and pk-2(y) denote a polynomial of degree 2k and (k — 2), respectively. We note 
that, as in (|3.10p . the integral of (u — uf) x against pk-2(y) is zero. Thus, surprisingly simple, 
we have 



(u - u I ) x ip x (x,y)dxdy 



K 



(u - tt/) x (4- + i i; l (y)^ X yk-^(x,0) + s^ 1 (y)ip xyk (x,0))dxdy 



(k+2). 



K 



1 r 



-1 



« - «/)z(4-i (y)^/=-i(^0) + s y ^''{y)i) xyk (x,^)) 



(*+!)/ 



y=l 
y=-l 



dx 



(u - u/) x . y (4_i(y)^fc-i( a; >0) + 4 {y)4' xv k{x^))dxdy. 



K 



(3.14) 



Let us consider the first boundary integral in (|3.14p . on the top edge of the square K. By (|3.4p 
and (133]) . 



(it - U/)a;(x, 1)4^1 WVW*- 1 ^ )^ 



(it - u/)(x, l)4_i (1)^-1 (s, 0) 



I=-l 



•ftp) 



(u — uj)(x, l)ip x 2 y k-i (x, 0)dx = 0, 



(3.15) 



noting again that tp x 2 y k-i (x, 0) is a Pfc_i polynomial in x only. The other boundary integral 
in (|3.14p is also as ip x 2 y k (x, 0) € Pfc-i too: 



(it - ui) x (x, l)s[ k+1 \l)ip xy k(x,0)dx 
(u - u/)(x, 1)4* +1) (l)Vv 0) 

s r>d " 



i 

s=-l 



(u — ui)(x, l)^ x 2 y k(x, 0)dx = 0. 



That is the boundary integrals in (|3.14p are all zero. We repeat the integration by parts in 
this direction, while the boundary terms would be zero by (|3. 12|) and (|3.5p . By the integration 
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by parts k times more, (|3. 14[) would be 

(u - ui) x il> x dxdy 



K 



K 



(u - ui) xy (s^\ip xy k-i (x, 0) + s^ +r> -ijj h (x, 0))dxdy 



(u - u I ) xy 2(s^_ ] L) i) xy k-i(x,0) + 3^i/) xy k(x,0))dxdy 



-i 



(ti-tij^Cx,!)^-! ^«fc-l(^0) + 4 ^xy"{x^))l = \dx 



_„ (n - u/) a , 3/2 (4 fc _ 1 1) V ; xy'=- 1 ( 2; 1 ) + s^'ip xy k(x,(}))dxdy 



(*). 



+ /„ ( u ~ u l)xy^ (sfc-iVV 2 - 1 (z, 0) + 0))dxdy. 



K 



xy 



->xy" 



(3.16) 



We will perform the integration by parts one last time. But this time, we will treat the two 
terms in the last integral differently. 



(u — ui) xy k+iSk-itp xv k-i(x,0)dxdy 



K 



xy" 



A 



(u — uj) xy k+is k ip xy k(x,0)dxdy 



K 



(u — ui) x 2 y k+iSk-iip y k-i (x, 0)dxdy 
^ [{u - u^^k+i Sk-itp y k-i (x, 0)] dy, 
(u — uj) xy k+2Skip xy k(x,0)dxdy. 



For the second integral, the boundary term disappears by the condition (|3.12p . For the first 
integral, we note that the boundary integrals will be cancelled due to the opposite line integrals 
on two sides of the vertical edge x = Xi or due to the boundary condition on ip. We note also 
that the (fe+ l)st and (fe + 2)nd partial derivatives on uj above are all zero. Hence, we get (|3.8p 
by summing over the estimation on all rectangles K € 7ft, plus a scaling and the Schwartz 
inequality, 



< 



(u- ui) x if> x dx. = V] / (u - u/^Vxdx 

K K 

l) fc+2 / (u x 2 y k+iSk-iip y k-i + u xy k+2Skip xy k) dx 

K 



K 



K 



<Ch k+2 \u\ H k+a {n) \ij\ H ^ n) . 

We note that the semi .f/^-norm is needed above to bound tp y k-i. Thus k > 1 is required. | 

In the proof, we can see that the decrease of one degree polynomial in y does not change 
the super-approximation of Qk+i,k in % direction. After (I3.16D . if we skip the last step of 
integration by parts, we would get the following corollary. That is, we avoid \ip y k-i \l 2 when 
k = 1 which cannot be bounded by \iP\h 1 - 
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Corollary 3.1 (one-order superconvergence) For any Qk+ik function ip € Vh,i, defined in 
(pTTI) . for any u € H k+2 {£1), and for all k>l, 

| / (u - ui) x tp x dx\ < Ch k+1 \\u\\ Hk+ 2U\\ H i. (3.17) 



Symmetrically, switching x and y in Lemma 13. 11 we prove the following lemma. 

Lemma 3.2 (two-order superconvergence) For any Qk,k+i function ip € Vh,2, defined in (|3.2p . 
and for any u G H k+3 (£l), if k > 1, 



(u- Ul ) y ^ y dx\ < Ch k+2 \\u\\ Hk+ 3\\^\\ H i. (3.18) 



For the same reasons in Corollary 13. 1\ we get the following corollary from Lemma 

Corollary 3.2 (one-order superconvergence) For any Qk,k+i function i/j € Vh,2, defined in 
(I3T2D . for any u G H k+2 (Q), and for all k > 1, 

| / (u-u I )yipydx\<Ch k+1 \\u\\ H k+2\\i}\\ H i. (3.19) 
Ju 



Though the interpolation order is (k + 2) in above two lemmas, only the (fc+1) order in two 
corollaries can be achieved in computation due to the coupling of terms in mixed formulation. 
We prove the approximation properties in the lower polynomial direction next. Now, even for 
k = 1, we have a two-order superconvergence. 



Lemma 3.3 (two-order superconvergence) For any Qk+i,k function ip € Vhi, defined in (13. ip . 
for any u € H k+3 (Q), and for all k > 1, 

(u - Ul ) y ip y dx\ < Ch k+2 \\u\\ Hk+3 \\iP\\ H x. (3.20) 

n 

Proof. Again, we first consider the estimation on the reference element K = [— 1, l] 2 . Since 
the polynomial degree in y is too low, we do Taylor expansion in x direction, different from 
the last lemma. 

ip y (x, y) = ^(0, y) + xif> xy (Q, y) + • • • + —ip xky (0, y) + rj; x k+i y (0, y). 

Again, similar to (|3.9p . the integral of {u — ui) y against x 3 terms are zero if j < k — 1, 
(u - ui) y x j ip x jy(Q, y)dxdy 



K 

l 



[(u - u^x^^yiO, y)] y x dx - l(u- ui)x^ x ,y2 (0, y)dxdy = 0, 



-1 " - y ~~ * JK 
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noting that x^tp x j y 2(0,y) € Qk-i,k-2- Using the polynomial function Sk(x) defined in (|3.1ip 
we have, cf. (I3.14p . 



(u — ui)yip y dxdy 



(u - n/)^(4 fe+2) (x)ip x k y {0,y) + (x)i/} x k+i y (0,y))dxdy 

' \„, / n \ , (fc+2) 



(k+3), 



-1 



(u-ui)y(s\ '(x)ip x k y (0,y) + s ( k+1 >{x)ip xk+ i y {Q,y)) 



x=l 



x=-l 



dy 



(fc+2). 



it- 
Here, for the first time integration by parts, the boundary integral disappeared by (|3.4p . 
(u — itj)(±l, ±1) = 0. In the next (k + 1) times of integration by parts, the boundary integrals 
on x = ±1 would be zero, directly by the boundary condition (|3.12p of Sfc(x). 

(u - u^yipydxdy = (-l) k+2 J (u - ur) x k+2y{s k ip xky (0, y) + s' k+1 ip x h+i y {0, y))dxdy. 

Thus, 



K 



(u — ui) y ipydxdy 



< Clinch 



< C\u\ 



+2 y\\L' 2 {k)\\' l i ) y\\L' 2 (k) 



iffc+3(X)l ( f ; ll^l(X)- 



The rest proof repeats that of Lemma 13.11 

As for above lemmas and corollaries, we can get the following corollary from Lemma 

Corollary 3.3 (two-order superconvergence) For any Qk+ik function ip E Vh \, defined in 
(fBTTj) . for any u € H k+2 {tt), and for all k > 1, 



(U — Ul)y1pydx\ < C/l fc+1 ||u|j i ^fe + 2||^||H 1 - 



(3.21) 



Corollary 3.4 For any Qkk+l function ip 6 Vh 2, defined in (|3.2p . and for any u £ H k+3 (Q,), 
and for all k > 1, 



/ (U — Ul)y1pydx\ < C/l fc+2 11^11^ + 3 || V" || H 1 , 

Jn 

(U — Uj)yTpydx\ < C/l fe + 1 ||li||^fc + 2||^'||^l. 



(3.22) 
(3.23) 



Now we study the superconvergence in the both bilinear forms. 

Lemma 3.4 For any (vh,Qh) £ x P^, defined in (|2.3p and (|2.4|) . and for any u S H 3 (Q,) (~1 
H l (fl), 

\a(u-I h u,v h )\ < Ch k+2 \\u\\ H k+ 3 ( n) 2\\v h \\ H i( n) 2, k>l, (3.24) 

\a(u - I h u, v h )\ < Ch k+1 \\u\\ Hk +2( Q) 2\\\r h \\ H i {n) 2, k>l, (3.25) 

1 6(u - I ft u, q h ) | < Ch k+1 \\u\\ Hk+ 2 {n) 2\\q h \\ L *(Q), k>l, (3.26) 
where I^u is the interpolation of u defined by (13. 3p . 
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Proof. (HEMP is a combination of (^g]) . ([3T20J) . (13T22D . and (13T8D . ([3T25D is a combination of 
§3TB), (IX2T]1 . (J3I23D, and (ETT9]) . 

For (|3.26f) . we will lose one order of convergence. Let qh = divw/j for some w/j = ((f>,ifi) G 
Vft. We have, denoting u = (u, f), 

6(u - I h u, q h ) = V" / ((« - tt/)a. + (w - vi) y )((j) x + ^ y )dx. 

Here we have two old integrals, /^-(m — uj) x <f> x dx and /^(f — vj) y ip y dx, and two new integrals, 
f K (u — uj) x ifj y dx and /^(f — vi) y (f) x dx. The approximation order can be one order higher for 
the two old integrals. For the two new integrals, by symmetry, we consider f K (u — ui) x Tp y dx.. 
We use the following Taylor expansion on the reference element K in the y direction. We 
note that the Taylor expansion in x direction would lead to a too high order polynomial in y 
direction each term in (|3.27p below. 

ip y (x,y) = ip y (x,0) +y?p y 2(x,0) H h ^ _ -y ip y k(x,0) + — ip y k+i(x, 0). (3.27) 

Here all ij} y j (x, 0) are polynomials of degree k in x. That is, a generic term y J tp y j+i (x, 0) G Qk,j- 
This is the same as the generic term y 3 ip xy j(x,0) in the early Taylor expansion (j3.9j) . Thus 
repeating the proof of Lemma 13. 1\ we get 

(u - ui) x ip y dx = / (u - ui) x (s^^ip k(x,0) + s k k+2 ' > ipk+i(x,0))dxdy 
k Jk 

= (-l) k+1 J u xy k+i(s k _^ y k(x,0) + s' k ip y k+i(x,0))dxdy. 

For the second integral, we can do an integration by parts to raise one more order. But we are 
limited by the first integral above to get only 



(u - U!) x ljj y dx 



< \\u\ 



Similarly, we have the same bound for \f^-(u — uj) y ^ x d^\. (13. 26|) follows by the Schwartz 
inequality and the scaling of referencing mappings. I 

Finally, we estimate the approximation to p. 

Lemma 3.5 For any function G ~Vh, defined in (|2.3p . and for any p G H k+1 (Q) n Lq(J1), 

I divv h (p-pi)dx\<Ch k+1 \\v h \\ H i\\p\\ Hk+ i, (3.28) 
Jn 

where pi is a special nodal interpolation of p in Ph, defined in (I3.29P below. 

Proof. We note that Ph are discontinuous Qk functions, Ph = div V^. We define an interpo- 
lation operator for Ph via that Ih for defined in (|3.3|) . For a p G H 2 (Q) n Lq(O), Arnold, 
Scott and Vogelius shown in [3] that there is a w G H 3 (Q) 2 n Hq(Q) 2 , such that 

divw = p, and ||w||#3 < C||p||#2. 
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For simplicity, we assume the above lifting exists up to order k + 1. We define 

pi = divwj, 

for w/ = I/jW defined by (|3.3j) . In order to use (|3.26p . we use notations: 



W=, ")' W/= W' Vfc \c 



Repeating the proof in Lemma 13.41 we get 



/ div Vh(p — pi)dx\ = | / divv/j div(w — w/)dx| 
Jn Jn 

((U - Ul) x + (V - Vl)y)((j) x + i>y)<ht 



in 

< Ch k+1 



H k+2 



1> 



H 1 



< Ch Upline hh\\m ■ 



(3.29) 



We derive the main theorem. 

Theorem 3.1 The finite element solution (\ih,Ph) of (|2.5p has the following superconvergence 
property, one order higher than the optimal order, 

\\ u h ~ Iftu|| H i + \\p h - pi\\ L 2 < Ch k+1 (\\w\\ H k+2 + \\p\\ H k+i), (3.30) 

where the interpolations (IhU,pi) are defined in (|3.3[) and (13.29H . 

Proof. By the inf-sup condition shown in |14} 125]. it follows that, cf. [T7], for all (w/^r/J € 
V,, x /',,. 



sup 

(v h ,q h )£V h xP h 



Chilli 1 + WML 2 

By Corollary 13.41 and Lemma 13.51 we have 



(Wfe, Vfc) + b(y h , r h ) + b(w ft , q h ) 

ii ij — |] — n > C(\\w h \\ H i + \\r h \\ L 2). (3.31, 



u h - IftuH^i + \\p h ~Pi\\l 2 

a(u h - IfeU, v fe ) + b(v h ,p h - pi) + b(u h - I h u, q h ) 



< C sup 

(v h ,q h )€V h xP h 



= C sup 

(v h ,q h )eV h xP h 

< Ch k+l 



\\ w h\\m + Wh\\L 2 

a(u-I h u,v h ) + b(v h ,p -pi) + b(u-I h u,q h ) 



IWhWm + Mh\\L 2 

\u\\ H k+2 + \\p\\ H k + l). 

Note that, due to the pointwise divergence free property, we have above that 

b(u h - I h u, q h ) = b(-l h u, q h ) = b(u - I h u, q h ). 
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Here, to be precise, we do not have a superconvergence for p in Theorem 13.11 As Ph 
are degree-/? polynomials, the best order approximation to p in L 2 -norm would be (k + 1). 
However, due to the mixed formulation, the convergence of ph to p is limited to the optimal 
order convergence of u^, which is (k — 1) in i? 1 -norm as has polynomial degree k only 
in y direction for its first component. In this sense, the superconvergence result (|3.30p does 
lift the order of ph by one. For k > 1, we may have two order superconvergence for the 
velocity. Such numerical examples are shown in [25] and in next section. That is, for some 
special functions u, I/jU might be also in the divergence- free subspace of V^. If so, we have a 
two-order superconvergence result. 

Theorem 3.2 (two-order superconvergence) For some solution u of (|2.ip . if 

I h u G Z h := {z h G V h | divz fc = 0} , 

where Ih is defined in (13.3H . and if k > 1, then 

\\u h - I h u\\ H i < C7t fc+2 ||u|| Hfc+ 3. (3.32) 

Proof. By (I3.24p . limited in to the divergence- free subspace, 

n T n a(u h -I h u,w h ) 
\\u h - Ihu\\ H i < sup 1| ii 

a(u-I h u,w h ) ^, fc +2|| I, 
= sup 1| < Ch T ||u||jjfe+3. 



4 Numerical tests 

In this section, we report some results of numerical experiments on the Qfc+i^-Qfc^+ielement 
for the stationary Stokes equations ()2. 1|) on the unit square Q, = [0, l] 2 . The grids Th are 
depicted in Figure El i.e., each squares are refined into 4 sub-squares each level. The initial 
grid, level one grid, is simply the unit square. 

We choose an exact solution for the Stokes equations (12. lj) : 

u = curl<7, p = Ag. (4.1) 

Here 

5 = 2V-z 4 )V-/) 2 . 

So we can compute the right hand side function f for (12. ID as 

f = -A curl # + VA#. (4.2) 

We note that, unlike [251 I14j . we intentionally choose a non-symmetric solution so that no 
ultraconvergence would happen, which does not exist in general. The solution p is plotted in 
Figure HI 
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Figure 4: The solution p (the errors are shown in Figure 0) 

We compute the Stokes solution on refined grids, cf. Figure El by the divergence Qk+i,k~ 
Qfc^+ielement (|2.3p and by the rotated Bernardi-Raugel element [5| [9l [T7]: 

Vh R = {vh € C(Q) 2 n H^V) 2 | v h \ K G Q k+lyk x Q fc>fe+1 VK € T ft } , (4.3) 
P^ R = {q h G Lg(n) | q h \ K G Q fe _i} . 

Following the analysis in p3], the stability of the rotated Bernardi-Raugel element would be 
proved. For the rotated Bernardi-Raugel element, the system of finite element equations is 
solved by the Uzawa iterative method, cf. [9l \T7\ I13j . The stop criterion is the difference 
\Ph^ ~ Ph ^1 — 10 -6 - We list the number of Uzawa iterations in the data tables by #Uz. 
Here the interpolation operators are standard Lagrange nodal interpolations 



Table 1: The errors = u — I^u and = p — pj for (I4.ip . 





\ e h L 2 


h n 




h n 


\\ e h L 2 


h n 






Qfc+i,fc-Qfc,fc+idivergence-free element, k = 1 


#it 


2 


0.264345 




1.341770 




5.965379 




4 


3 


0.102329 


1.4 


0.795594 


0.8 


1.896372 


1.7 


4 


4 


0.026839 


1.9 


0.219469 


1.9 


0.481076 


2.0 


3 


5 


0.006773 


2.0 


0.055901 


2.0 


0.120363 


2.0 


3 


6 


0.001697 


2.0 


0.014035 


2.0 


0.030083 


2.0 


3 


7 


0.000424 


2.0 


0.003512 


2.0 


0.007520 


2.0 


3 




rotated Bernardi-Raugel element (|4.3j) 


#Uz 


2 


0.570990 




3.531380 




7.497615 




29 


3 


0.244967 


1.2 


3.028368 


0.2 


6.943183 


0.1 


65 


4 


0.074335 


1.7 


1.797533 


0.8 


3.300598 


1.1 


136 


5 


0.019849 


1.9 


0.946426 


0.9 


1.575390 


1.1 


297 


6 


0.005080 


2.0 


0.481087 


1.0 


0.762341 


1.0 


330 


7 


0.001281 


2.0 


0.241916 


1.0 


0.373990 


1.0 


204 



For the (5fc+i,A;-Qfc,A;+idivergence-free element, the pressure does not enter into computation, 
but is obtained as a byproduct, because Ph = divV^. The resulting linear system of Qk+i,k- 
Qfc^+idivergence-free element equations can be formulated as symmetric positive definite. 
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Then the iterated penalty method |13[ [25] can be applied to obtain the divergence-free finite 
element solution for the velocity, and a byproduct ph = divw/j for the pressure. In our 
computation, the iterated penalty parameter is 2000. The stop criterion is the divergence 
|| divu^||o < 10~ 9 . The number of iterated penalty iterations is also listed as #it in the data 
tables. 

In Table [H we list the errors in various norms for the (5fc+i 5 ArQfc,fc+idivergence-free element 
and for the rotated Bernardi-Raugel element, for k = 1. It is clear that the order of convergence 
is 2, one order higher than that of latter. We note that the convergence order is only 2 for 
Q2,i-Qi,2 divergence-free elements in L 2 -norm, i.e., no L 2 -superconvergence. But we do see 
L 2 -superconvergence for k > 1 next. 



Table 2: The errors = u — I/jU and = p — pj for (|4. 1 j) . 





\ e h L 2 


h n 




h n 


l e /ilU 2 


h n 






Qfc+i,fc-Qfc,fc+idivergence-free element, k = 2 


#it 


1 


0.322530 


0.0 


1.580066 


0.0 


3.546405 


0.0 


3 


2 


0.071851 


2.2 


0.699614 


1.2 


1.010498 


1.8 


4 


3 


0.005510 


3.7 


0.089611 


3.0 


0.131816 


2.9 


3 


4 


0.000355 


4.0 


0.010471 


3.1 


0.015587 


3.1 


3 


5 


0.000022 


4.0 


0.001280 


3.0 


0.001904 


3.0 


3 


6 


0.000001 


4.0 


0.000159 


3.0 


0.000236 


3.0 


3 


7 


0.000000 


4.0 


0.000020 


3.0 


0.000029 


3.0 


3 




rotated Bernardi-Raugel element (j4.3j) 


#Uz 


1 


0.645475 


0.0 


4.250791 


0.0 


1.143688 


0.0 


27 


2 


0.191342 


1.8 


2.518701 


0.8 


5.499136 




67 


3 


0.025892 


2.9 


0.673622 


1.9 


1.621194 


1.8 


100 


4 


0.003307 


3.0 


0.172036 


2.0 


0.441596 


1.9 


156 


5 


0.000419 


3.0 


0.043543 


2.0 


0.113029 


2.0 


266 


6 


0.000053 


3.0 


0.010954 


2.0 


0.028424 


2.0 


130 


7 


0.000007 


3.0 


0.002747 


2.0 


0.007117 


2.0 


101 



In Table [21 we list the computation results for k = 2 elements. Again, the divergence-free 
element is one order higher than the rotated Bernardi-Raugel element. To show the difference 
in the two elements, we plot the errors by two elements on level 4 grid in Figure [5l One 
can see the advantage of the divergence- free element, which fully utilizes the approximation 
power of Ufr by lifting the pressure polynomial degree. Of course, another advantage is the 
divergence-free solution after such a lift. We finally report the results for k = 3 in Table [3j 
All numerical results confirm the theory, and also show the sharpness of the super convergence 
analysis. 

Finally, we test the two-order super convergence in Theorem 13.21 We choose a symmetric 
function as the exact solution of the Stokes equations (|2. 1 j) : 

u = curl 5 , g = 2 8 {x-x 2 ) 2 {y-y 2 ) 2 . (4.4) 

Comparing to the data in Table El we can see, in Table HI that the velocity does converge with 
another order higher than the optimal order. This is predicted in (|3.32p . Here the order of 
convergence for the pressure is the same as that in Table [3l It indicates that the analysis in 
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Figure 5: The errors of ph for the divergence- free (top) and BR elements. 

Theorem 13.11 is sharp. Here we have an order-two super convergence in L 2 -norm too, for the 
velocity. But this is not proved in this manuscript. 
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| e /i |l 2 




\ e h\m 
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\ e h\m 
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